function element = computeRenh(e,element)
crackLength = element.crackLength;
J = computeJ(element);
u_mu = J*element.Xi;
u_mu_n = u_mu'*element.normal;
T = element.ft + element.H*abs(u_mu_n);
Tn = max(T,0);
Tm = 0;
TGamma = Tn*element.normal + Tm*element.tangent;

%% signalize if element's crack got active
if T<0
    element.active = 0;
else
    element.active = 1;
end

%% compute renh
t = element.thickness;
renh = -crackLength*J'*TGamma*t;
ng = e.getNumGaussPoint();
f = e.getFace();
for j=1:ng
    xi = e.getGaussPoint(j-1);
    w = e.getGaussWeight(j-1);
    jac = f.jacXAt(xi);
    detJac = det(jac);
    Ge = computeGe(element,xi);
    renh = renh - Ge'*...
        element.stress(:,j)*(w*detJac*t);
end
element.renh = renh;

%% update Kxixi if element change crack state
if element.active == 0 && element.active_old == 1
    %fully opening crack/crack close
    element.Kxixi = element.Kxixi_1;
    element.active_old = 0;
    disp(['***warning: crack is not active at element ' num2str(f.getId())]);
elseif element.active == 1 && element.active_old == 0
    %partial opening crack
    Ctilde = sign(u_mu_n)*computeCtilde(element);
    element.Kxixi = element.Kxixi_1 + ...
        crackLength*J'*Ctilde*J*t;
    element.active_old = 1;
    disp(['***warning: element ' num2str(f.getId()) ' got active crack']);
end
